Mesoscopic Model for Free Energy Landscape Analysis of DNA sequences. 
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A mesoscopic model which allows us to identify and quantify the strength of binding sites in DNA 
sequences is proposed. The model is based on the Peyrard-Bishop-Dauxois model for the DNA chain 
coupled to a Brownian particle which explores the sequence interacting more importantly with open 
base pairs of the DNA chain. We apply the model to promoter sequences of different organisms. 
The free energy landscape obtained for these promoters shows a complex structure that is strongly 
connected to their biological behavior. The analysis method used is able to quantify free energy 
differences of sites within genome sequences. 



PACS numbers: 87.15.A-, 87.15.H-, 87.14.gk 

I. INTRODUCTION 

The study of biomolecules at a mesoscopic level tries 
to identify the main degrees of freedom of the system and 
understand its behavior in terms of the dynamical and 
statistical-mechanics properties of the model. At this 
level, the concept of a free-energy landscape (FEL) rep- 
resents a paradigm for the comprehension of several bio- 
logical complex problems such as protein folding, protein 
structure, and biomolecular interaction [1 j. A FEL gives 
the change in the free energy of a system when the dif- 
ferent degrees of freedom change. The description given 
by the topology of the FEL permits to connect struc- 
ture, dynamics, and thermodynamics in many different 
systems ranging from atomic clusters to biomolecules or 
soft matter systems [2]. 

Here we address our attention to the characterization 
of the FEL of DNA sequences and, in particular, those 
sites with regulatory and transcriptional relevance. Re- 
cently, a great attention has been devoted to the mecha- 
nism whereby proteins bind to specific sites on DNA [3]. 
The quantification, grounded in a physical basis, of the 
strength of these binding sites is an open problem. In this 
paper, we propose a model which allows us to calculate 
free energy differences between specific and nonspecific 
binding sites. Even more, we are able to build, from 
the trajectories of the model, a representation of the free 
energy landscape. 

The model is inspired in the protein search of the bind- 
ing sites in a DNA chain. This search is a combination 
of three-dimensional (3D) jumps between separated re- 
gions of DNA and one-dimensional (ID) diffusion along 
the chain. It has been stated that most of the search 
time is spent in the ID diffusion process since the time 
jumps in three dimensions is negligible |3HB]. Thus, the 
restriction to a ID search is a good starting point for our 
model. On the other hand, it has also been conjectured 
that the dynamics of the DNA chain plays an important 
role in the recognition of binding sites by the regulatory 



factors or the transcription protein [7J [8]. Thus, tran- 
scription processes, for instance, would be induced by 
the binding of RNA-polymerase to openings (bubbles) 
in the DNA chain. This idea is supported by compu- 
tational approaches to the DNA dynamics [9 , and ex- 
perimental evidences [10 j. Our model follows this idea 
and considers the interaction of a test particle, which ex- 
plores the DNA chain, coupled to the bubbles. In order 
to be physically and biologically relevant, such a descrip- 
tion should provide useful qualitative and quantitative 
information about the process. A similar strategy has 
been used to characterize complex networks and identify 
regions of special relevance (communities). In such an 
approach, a "fictitious" Brownian particle goes over the 
graph [11] and its dynamics reveals "thermodynamics" 
and structural quantities of the topology [T2] . 

In this paper, we propose a mesoscopic model for DNA- 
particle interaction. In our picture, the test particle un- 
dergoes a ID Brownian motion in interaction with a clas- 
sical field, the DNA chain itself, whose dynamics is also 
affected by the presence of the particle. The test particle 
interacts more strongly with open base-pairs of the DNA 
chain. In this way "softer" regions of the DNA sequence 
are more likely to be visited by the particle, which will 
help also in stabilizing the bubbles. This interaction is 
not sequence-dependent, as the DNA base-pair dynam- 
ics already depends on the sequence. Thus, this model 
could also represent the interaction of a real protein as 
RNA-polymerase, with the DNA bubbles. 

Particle and chain are described at the same level of 
complexity. We use the Peyrard-Bishop-Dauxois (PBD) 
[T3j [14] model to perform the dynamics of the chain. This 
model was proposed initially for the study of DNA ther- 
mal denaturation and incorporates the formation and dy- 
namics of bubbles in a natural way. The PBD model can 
take into account the sequence information through its 
parameters. 

Our model incorporates three basic ingredients of the 
physics of the system: DNA sequence, bubble dynam- 
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FIG. 1. (Color online) Schematic illustration of the DNA- 
particle interaction model. The unidimensional chain (solid 
circles) is represented by the base pair opening coordinate (yi) 
while the probe particle (shaded ellipse) is a diffusing particle 
along the DNA chain (X p ). 



ics, and ID particle diffusion. The analysis of the DNA- 
particle complex allows us not only to identify possible 
binding sites, but also to describe the whole structure of 
the free-energy landscape and determine free-energy dif- 
ferences between different representative states. We show 
the validity and usefulness of our approach by studying 
the FEL of three promoter sequences. In each case, the 
FEL topology gives insight into the biological behavior 
of the system. 



II. MODEL 

We describe the DNA chain by a modified Peyrard- 
Bishop-Dauxois model [T3l - [T5] . There, the complexity of 
the DNA molecule is reduced to the study of the dy- 
namics of N base pairs described by the variables y n 
defining the distance between the bases. In this frame- 
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where the first term accounts for the hydrogen bond in- 
teraction and the second one for interactions with the 
solvent [15] [16]. It was shown in [15 that the inclusion 
of this barrier modifies drastically the duration and sta- 
bility of the bubbles. 

The potential W(y n ,y n -i) describes the stacking in- 
teraction between the base pairs along the DNA strand, 



W{y n ,y n - 1 ) = ±K(l + pe- s to« 
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In order to study different DNA sequences, the PBD 
model includes sequence dependent Morse parameters: 
D ni a n . Regarding the DNA chain, we use the set of 
parameters [17] considered in [15] . 

The particle is represented by a Brownian particle (see 
Fig. [T]) moving in a one-dimensional space with coordi- 
nate X p and interacting with the DNA chain through a 
phenomenological potential which depends on X p and the 
set of coordinates {yi}f =1 : H P = pl/2m p + V int (X p , {yi}) 
with 
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where B sets the interaction amplitude, 7 the range of 
interaction with the base separation and a the spatial 
range of interaction on the DNA chain. The functional 
form for the interaction has been chosen to be linear at 
low yi and to saturate at large yi in order to avoid that 
the chain opens indefinitely. Note that with this term the 
particle is trying to open the chain in a length range of a 
and get self-trapped. Although possible, no sequence de- 
pendence is included in this term since we are interested 
in giving a general picture of the FEL. 

We still have to fix the parameters for this interaction 
term. For the particle damping and mass, we take r\ v 
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and mp = 7000Da. These values are of the 



order of magnitude of proteins which bind DNA [T8l [19] . 
The intensity of the interaction chosen is B = 0.52eV. 
This value provides local interactions of the order of the 
Morse potential dissociation energy at each base pair. 
The parameter 7 = 0.8 A -1 saturates the interaction at 
y = 1.25 A, a typical value for open base pairs. We take 
a = 1 for the longitudinal separation between base pairs, 
in arbitrary units, and consider a = 3. This provides 
an interaction range of around 5 or 6 base pairs (bp). 
It is interesting to note that this value has been chosen 
in order to observe states with bubbles of 10 — 20 bp, 
which is an adequate width for the processes we take 
into account here [20] . 

Once we have fixed the model parameters, we derive 
the Langevin equations for both the chain bases and the 
particle. For the chain, we get: 
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where r\ stands for the damping and £ ra for the thermal 
noise, so (£„(£)) = and = 2rm]k B T6 nk 5(t- 

f) hold. 

The Langevin equation for the particle is: 



d 2 X p 
dt 2 



m p r]p- 



dV ir 



dt 



dX r , 



+ &(*)■ (5) 



where, r] p stands for the particle damping and £ p for the 
thermal noise. The fluctuation-dissipation relation reads 
as (£ p (t)) = and = 2m p r] p k B TS(t - t'). 

The equations were numerically integrated using the 
stochastic Runge-Kutta algorithm [21 . The integration 
of the Langevin equations of motion provides trajectories 
of the particle and the DNA chain. Each DNA sequence 
was simulated in five different realizations for 40/is, us- 
ing 10 fs time steps and a lfis preheating time. Since 
it has been reported that ID diffusion periods cover a 
time of the order of milliseconds [18] [19], the simula- 
tion time used is reasonable for the problem considered 
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FIG. 2. (Color online) Trajectory of the particle shown in red 
line (light gray) with the DNA chain for a collagen sequence 
(see results section). Black points stands for open bp (vat > 
l.A or ycG > 1.5 A) and white for closed ones. 



here. The simulation temperature is T = 290K and the 
boundary conditions for the protein are periodic, while 
for the chain we consider the hard boundary conditions 
discussed in [15 . An example of such a trajectory is 
given in Fig. |2j It is observed that the particle moves 
in a "sea" of open bubbles which clearly shows the soft 
domains of the genome structure. The particle eventu- 
ally jumps between these domains opening large stable 
bubbles. As seen in the figure the dynamics of the bases 
is strongly affected by the presence of the particle. 



III. ANALYSIS 

To extract useful information of such trajectories, and 
due to the large dimensionality of the system, we ap- 
ply the principal component analysis (PCA) [22] to the 
chain trajectory. It has been proved [15 that the first 
few eigenvectors unveil the softest regions of the DNA 
chain and hence the possible binding sites for our parti- 
cle. Even more, PCA reduces the large number of degrees 
of freedom of the system to just a few, by projecting the 
coordinates of the system into the first few eigenspaces 
(reduced trajectories). For each of the sequences consid- 
ered, we restrict ourselves to the first five eigenspaces. 
This subspace accounts for 75% of the total fluctuations 
of the chain dynamics. 

To obtain the FEL properties of the system we make 
use of the map of trajectories to a conformational Markov 
network (CMN) [17,20]. The CMN has been proven 
to be a useful representation of large stochastic trajec- 
tories [24-26 . This coarse grained picture is usually 



constructed by discretizing the conformational space ex- 
plored by the dynamical system and considering the hops 
between the different configurations as dictated by the 
MD simulation. In this way, the nodes of a CMN are 
the subsets of configurations defined by the conforma- 
tional space discretization, and the links between nodes 
account for the observed transitions between them. The 
information of the stochastic trajectory allows us to as- 
sign probabilities for the occupation of a node (Pi) and 
for the transitions between two different configurations 
(Pij). Defined as above, a CMN is thus a weighted and 
directed graph. It should be stressed that the informa- 
tion contained in the CMN is much richer that one given 
by equilibrium statistical mechanics since it includes the 
dynamics of the system encoded in the probability tran- 
sitions, P^. 

In our case, we start from the reduced trajectory for 
the DNA (obtained using five principal components) and 
the trajectory of the test particle. We discretize the to- 
tal coordinate space in 20 bins of equal volume for the 
reduced trajectory and N bins (the DNA base pairs) for 
the particle. This will constitute the microstate space 
of the CMN, each node with occupancy probability Pi 
obtained from the reduced trajectory. Once the CMN 
has been built, we split it into basins of attraction, i.e., 
regions in which the probability fluxes (Pij) converge to 
a common state (attractor) of the network. This task is 
usually hard, since algorithms scale as power law of the 
system size. In this case we have applied the stochastic 
steepest descent algorithm developed in [23] , which scale 
as N log N. In this decomposition, a basin corresponds to 
a coarse-grained state (of connected nodes) of the CMN. 
In next section, we represent each basin by its attracting 
node. 

Once these basins have been defined we can represent 
the FEL by a hierarchical tree diagram (dendrogram) 
[23] . built according to the weights and links among the 
basins. This representation is similar to the "disconnec- 
tivity graph" scheme used in other context [TJ [2j [27] . 
First, an "adimensional free energy" is assigned to each 
node i given by Fi/kT = log(P w ) — log (Pi) where w rep- 
resents the weightiest node. Using this magnitude as con- 
trol parameter, we slowly increase it step by step from 
its zero initial value. At each step of this process, we 
obtain a network composed of those nodes with free en- 
ergy lower than the current threshold value. As the free- 
energy threshold increases, new nodes emerge together 
with their links. These new nodes may be attached to 
any of the nodes already present in the network or they 
can emerge as a disconnected component. At a certain 
value of F I kT, some components of the network become 
connected by the links of a new node incorporated at 
this step. Initially we have a set of disconnected vertical 
lines (corresponding to basins) which become linked once 
the control parameter has overcome the barriers between 
them i. e. when the free energy of the saddle nodes is 
reached. Then we draw a horizontal line linking these 
two basins. Obviously, for large threshold all the net- 
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FIG. 3. (Color online) Free energy dendrograms (top) for each of the three sequences together with selected states (bottom). 
The dendrograms represented are a detail of the whole structure showing low energy basins only. The construction of the 
dendrogram has been performed considering the representative node of each basin. Significant biological states have been 
searched within the network structure. They correspond to the important basins of the hierarchical free energy organization. 



work is connected. We can plot this process as a "tree 
diagram" or dendrogram. 

Using this representation we can understand qualita- 
tively and quantitatively the hierarchical organization of 
the basins and the barriers among them and figure out 
the behavior of different sequences. 

This method could be applied to a DNA chain without 
a particle. However, the inclusion of the particle is es- 
sential to get the FEL of the system. An analysis of the 
DNA alone (PBD model) lets us determine the opening 
probabilities and average position of the chain base pairs, 
and unveils the softer regions that can indeed be related 
to sites of biological importance. Nevertheless, the FEL 
of this model is trivial, as opening events are rare and 
the chain remains closed for most of the time. The in- 
clusion of the particle stabilizes the bubbles (as can be 
observed in Fig. |2| and allows us to go further in terms 
of predictions. We are able to define relevant states in a 
precise and systematic way (basins), to predict possible 
binding sites, and to extract the thermodynamical mag- 
nitudes related with them, thus characterizing these sites 
in terms of biological importance. 

IV. RESULTS 

To illustrate the method and validate our model, we 
analyze three different promoter sequences. Promoters 
are DNA regions in which regulation and initiation of 
transcription of a gene occurs. Two of them correspond 
to the so called strong promoters, while the one left is a 
weak promoter [28]. Strong promoters show a high level 
of expression in mRNA and usually their sequences are 
close to the consensus sequence. The strong promoters 
studied here are the P5 virus promoter, given by the 69 



bp sequence shown in [7] and the human collagen type 
I a2 chain, given by the 80 bp sequence shown in [9]. 
Finally, the weak promoter is the lac operon regulatory 
region, whose 129 bp sequence has been taken from [8]. 

In Fig. [3] (top) we show a detail of the free energv 
dendrograms for each of the mentioned sequences. The 
basin structure consists of a big set of low occupied (high 
energy) basins, and a small set which gathers almost the 
whole trajectory (see below). This small set of basins is 
the one shown in Fig. [3j 

In the botton of Fig. [3] some remarkable states for each 
of the three promoters are highlighted. The method iden- 
tifies states with a biological meaning as they correspond 
to the most important basins. The most significant sites 
we are dealing with are the transcription starting site 
(TSS) and the TATA box, although additional promoter 
sequences can be found depending on the genome. The 
RNA polymerase binds to the TSS, starting the tran- 
scription into mRNA. Promoter sequences are usually 
labeled from the TSS (+1). The TATA-box is found ap- 
proximately 35 base pairs upstream from the TSS [28] . 

For the collagen chain, three states have been high- 
lighted. State A identifies the TSS, showing a bubble in 
this region with the particle placed just there. States B 
and C are linked to excitations of other important sites 
such as the TATA box (state B), see [9 . In the same 
way, we have found a basin related to the TSS in the 
case of the P5 chain (state C) and the lac operon (state 
C), together with other regulatory sites. 

The arrangement of these basins in the free energy den- 
drogram informs about the relative free energy between 
the states and the relation between them. For example, 
the collagen dendrogram contains three main branches, 
each one related to each of the three states shown. The 
P5 promoter shows an analogous structure: two main 
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branches and another one divided into two states (B and 
C) which are kinetically close. The remaining states of 
each branch correspond to states similar to that shown, 
with only slight variations in the chain conformation or 
in the particle position. 

The lac promoter shows a clearly different behavior 
compared with the two strong promoters. From a qual- 
itative point of view, the arrangement of basins differs 
from the P5 sequence or the collagen one. To visual- 
ize quantitatively the difference, the basin occupancy is 
plotted in Fig. [4j We show the weight of each basin 
(red bars) for the three sequences together with the ac- 
cumulated weight (blue line). It is remarkable, in the 
case of the collagen sequence, that a few basins (25 out 
of 1661) accumulate almost the whole weight of the net- 
work (over the 99%). The results for the P5 promoter 
are completely analogous, a few basins account for most 
information of the dynamics. These basins are the ones 
shown in the dendrogram of Fig. [3J When we inspect 
the bottom graph in Fig. [4] we see a completely different 
tendency. In the lac network, the distribution of weight 
among the basins is more uniform. 

For the collagen and P5 networks we can define a 
threshold from which the individual contribution to the 
total weight is negligible. The trajectory is concentrated 
in around 25 basins and the remainder of the network can 
be seen as a "background" . This "background" is limited 
to those basins with a weight below 10 -3 (horizontal lines 
in Fig [4]). Following this criterion, we can distinguish be- 
tween specific and nonspecific states. Those basins above 
the threshold (vertical lines in Fig. [4]) may be defined 
as specific states (with a clear biological function) while 
those below the threshold may be defined as non-specific 
states. 

For the collagen sequence, 25 "specific" basins appear, 
covering 99.41% of the total trajectory. The P5 and lac 
sequences show respectively 23 and 88 "specific" basins, 
which gather 99.38% and 96.91% respectively of the total 
network weight. Using these definitions, we are able to 
calculate the relative free energy between states. These 
magnitudes reveal the "strength" of the different sites in 
each promoter. It has been reported that specific bind- 
ing proteins show a greater affinity for strong promoters 
than for weak ones [29]. To quantify these differences 
we calculate thermodynamical properties of the most im- 
portant basins. Once we have divided the network into 
the different basins of attraction, several statistical mag- 
nitudes can be defined from them. The weight of the 
basin is defined as the sum of that of the nodes belong- 
ing to the basin, i. e. for a basin a we have P a = Pi 
with i G a. In the same way the entropy of each basin 
can be defined as S a /kB = — X^^UogPi with i G a. 
Attending to the previous definition of the nonspecific 
basin, whose thermodynamical magnitudes can be com- 
puted as explained, we can calculate the free energy of 
each basin with respect to the nonspecific state. If Pp 
is the weight of the non-specific basin, then the free en- 
ergy difference between a basin a and the macrostate 
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FIG. 4. (Color online) Basin occupancy (bars) together with 
accumulated weight (solid line) for the collagen, P5 and lac 
sequences. The horizontal line shows the weight threshold be- 
tween specific and nonspecific states. The vertical line estab- 
lishes such a frontier in terms of basins. Note the logarithmic 
y axis. 



P is AF a /k B T = -\og(P a /Pp). Table [!] shows signifi- 
cant differences between strong and weak promoters. On 
the one hand we observe that both the total weight and 
entropy of the non-specific states in the weak promoter 
exceed by almost an order of magnitude the ones shown 
for the strong promoters. On the other hand, we can see 
that the specific states show much higher free energy dif- 
ferences with respect to the nonspecific states in the case 
of the strong promoters than the ones shown for the lac 
sequence. Thus, the analysis presented here opens the 
way to a systematic study of promoter character within 
the framework of a mesoscopic model. 



6 



Promoter 


State 


Pe 


S/k B 


—AF/kT 




A (TSS) 


0.169 


1.365 


3.305 




B(TATA) 


0.157 


1.380 


3.232 


Collagen 


C 


0.086 


0.652 


2.519 




NS 


0.006 


0.085 


0.000 




A (TATA) 


0.135 


1.051 


3.130 




B 


0.107 


0.913 


2.898 


P5 


C (TSS) 


0.086 


0.684 


2.681 




D 


0.059 


0.494 


2.301 




NS 


0.006 


0.027 


0.000 




A (TATA) 


0.115 


0.970 


1.311 




B 


0.095 


0.891 


1.120 


lac 


C (TSS) 


0.090 


0.775 


1.066 




D 


0.038 


0.373 


0.204 




NS 


0.031 


0.390 


0.000 



TABLE I. Statistical (occupation probabilities) and thermo- 
dynamical (entropy and free energy differences) magnitudes 
calculated for the chosen states of Fig. [3] and the nonspecific 
states (NS), according to the criteria shown in Fig. [4] 



In addition to the three promoter sequences of real bi- 
ological systems, we have analyzed a random sequence 
in order to prove the validity of our model. The ran- 
dom sequence has been obtained by taking the P5 pro- 
moter sequence and shuffling its base pairs, so that the 
obtained sequence contains the same base pairs but in 
random positions. This sequence should contain no ge- 
netic information at all, and this fact must be reflected 
in our analysis. 

When analyzing the random sequence with our 
method, we can see huge differences compared with the 
P5 promoter, as we would expect (see Fig. |5|. First the 
structure of the network is completely different. As there 
are no prominent states in the sequence, the number of 
basins is huge (8388 compared with the 529 in the P5 
promoter). The first two eigenvectors are representative 
of a homogenous lattice without localized states. The 
distribution of weights is also clearly different as now 
the "background" basins keep 6% of the total network 
weight, an even higher value than that of the background 
basins gathered in the weak promoter. The dendrogram 
also shows a much more distributed structure where, even 
though some nodes appear to fall to privileged positions, 
their relevance within the whole network structure is far 
from being comparable to that shown in networks from 
biological promoters. All these facts validate our model, 
as we can clearly distinguish between a sequence with 
binding sites, and thus with biological information, and 
one with none, even though their chemical composition 
is the same. 



V. CONCLUSIONS 

In this paper we have proposed and analyzed a meso- 
scopic model for the characterization of binding sites on 
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FIG. 5. (Color online) Analysis of a random sequence. Top: 
Probability of aperture along with the first two PCA eigen- 
vectors. Middle: Free energy dendrogram. Bottom: Basin oc- 
cupancy (bars) together with accumulated weight (solid line) 



DNA promoter sequences. The model is based on the ID 
diffusion of an extended probe particle along the DNA 
chain. The particle is coupled to the opening states of 
the chain (bubbles). In its dynamics, it visits the main 
sites of the sequences, with dwelling times covering a high 
percentage of the trajectory. Such behavior has allowed 
us to perform a deep analysis of the FEL which reveals 
the structure of the complex phase space. The analyzed 
promoter sequences have been chosen to include genomes 
from organisms of different domains (virus, bacteria and 
eukaryote) and different strengths of expression. The 
model and the analysis used are able to capture the main 
biological details of the sequences. 

Our model gives energy differences between specific 
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and nonspecific sites of the promoter. Our results are in 
good relative agreement with some data in the literature 
(see for instance [29]): they account for energy ratios be- 
tween weak and strong promoters. This fact would also 
make possible the study of sequences in which several 
TSSs are involved, showing the relative strength between 
them. 

We think that our results show the power of coarse- 
grained or phenomenological mesoscopic models to qual- 
itatively and quantitatively analyze complex biological 
systems, in particular the problem of protein-DNA reg- 
ulatory and transcriptional interactions. Protein-DNA 
interaction is a fundamental problem which has been 
the object of a very intense research from many differ- 
ent points of view in the past years [3j [5] . Our system 



can be seen as the searching problem of a universal pro- 
tein on a given DNA sequence, providing an approach 
for the study of specific protein-DNA interactions at the 
mesoscopic level, where different protein will interact in 
different ways with DNA molecules. 
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